#include <cstdlib>
#include <cstring>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>

#include <itkImage.h>
#include <itkImageFileWriter.h>
#include <itkAddImageFilter.h>

const unsigned int Dim = 2;
typedef double TPixel;
typedef itk::Image< TPixel, Dim >                     TImage;
typedef itk::AddImageFilter< TImage, TImage, TImage > TAddFilter;
typedef itk::ImageFileWriter< TImage >                TImageWriter;

int main( int argc, char* argv[] )
{
  std::string input_prs_file = argv[ 1 ];
  std::string output_image_file = argv[ 2 ];

  std::ifstream input_prs( input_prs_file.c_str( ) );
  if( !input_prs )
  {
    std::cerr << "Error with " << input_prs_file << std::endl;
    return( 1 );

  } // fi

  char buffer[ 5000 ];

  input_prs.getline( buffer, 5000 );
  std::string str = buffer;
  if( str.find( "[M]" ) == std::string::npos )
  {
    std::cerr << "Error header of " << input_prs_file << std::endl;
    return( 1 );

  } // fi

  input_prs.getline( buffer, 5000 );
  str = buffer;

  input_prs.getline( buffer, 5000 );
  str = buffer;
  while( str.find( "[...]" ) == std::string::npos )
  {
    input_prs.getline( buffer, 5000 );
    str = buffer;

  } // elihw

  std::string strnCols = str.substr( str.find_last_of( "[" ) );
  strnCols = strnCols.substr( 1, strnCols.find_last_of( "]" ) - 1 );
  unsigned int nCols = std::atoi( strnCols.c_str( ) );

  input_prs.getline( buffer, 5000 );
  str = buffer;
  std::vector< std::vector< double > > pseudo_image;
  while( str.find( "[M]" ) == std::string::npos )
  {
    std::string line = str.substr( str.find_first_of( "]" ) + 1 );
    pseudo_image.push_back( std::vector< double >( ) );

    char tokenize[ 5000 ];
    std::strcpy( tokenize, line.c_str( ) );
    char* tok = std::strtok( tokenize, "\t\n " );
    while( tok != NULL )
    {
      double val = std::atof( tok );
      pseudo_image.rbegin( )->push_back( val );

      tok = std::strtok( NULL, "\t\n " );

    } // elihw

    input_prs.getline( buffer, 5000 );
    str = buffer;

  } // elihw

  TImage::SizeType size;
  size[ 0 ] = pseudo_image.size( );
  size[ 1 ] = nCols;

  TImage::SpacingType spacing;
  spacing.Fill( 1 );

  TImage::PointType origin;
  origin.Fill( 0 );

  TImage::Pointer img1 = TImage::New( );
  img1->SetRegions( size );
  img1->SetSpacing( spacing );
  img1->SetOrigin( origin );
  img1->Allocate( );

  std::vector< std::vector< double > >::const_iterator pIt =
    pseudo_image.begin( );
  for( unsigned int col = 0; pIt != pseudo_image.end( ); pIt++, col++ )
  {
    std::vector< double>::const_iterator rIt = pIt->begin( );
    for( unsigned int row = 0; rIt != pIt->end( ); rIt++, row++ )
    {
      TImage::IndexType idx;
      idx[ 0 ] = row;
      idx[ 1 ] = col;

      img1->SetPixel( idx, *rIt );

    } // rof

  } // rof

  input_prs.getline( buffer, 5000 );
  str = buffer;
  while( str.find( "[...]" ) == std::string::npos )
  {
    input_prs.getline( buffer, 5000 );
    str = buffer;

  } // elihw

  strnCols = str.substr( str.find_last_of( "[" ) );
  strnCols = strnCols.substr( 1, strnCols.find_last_of( "]" ) - 1 );
  nCols = std::atoi( strnCols.c_str( ) );

  input_prs.getline( buffer, 5000 );
  str = buffer;
  pseudo_image.clear( );
  while( !input_prs.eof( ) )
  {
    std::string line = str.substr( str.find_first_of( "]" ) + 1 );
    pseudo_image.push_back( std::vector< double >( ) );

    char tokenize[ 5000 ];
    std::strcpy( tokenize, line.c_str( ) );
    char* tok = std::strtok( tokenize, "\t\n " );
    while( tok != NULL )
    {
      double val = std::atof( tok );
      pseudo_image.rbegin( )->push_back( val );

      tok = std::strtok( NULL, "\t\n " );

    } // elihw

    input_prs.getline( buffer, 5000 );
    str = buffer;

  } // elihw

  size[ 0 ] = pseudo_image.size( );
  size[ 1 ] = nCols;

  spacing.Fill( 1 );

  origin.Fill( 0 );

  TImage::Pointer img2 = TImage::New( );
  img2->SetRegions( size );
  img2->SetSpacing( spacing );
  img2->SetOrigin( origin );
  img2->Allocate( );


  pIt = pseudo_image.begin( );
  for( unsigned int col = 0; pIt != pseudo_image.end( ); pIt++, col++ )
  {
    std::vector< double>::const_iterator rIt = pIt->begin( );
    for( unsigned int row = 0; rIt != pIt->end( ); rIt++, row++ )
    {
      TImage::IndexType idx;
      idx[ 0 ] = row;
      idx[ 1 ] = col;

      img2->SetPixel( idx, *rIt );

    } // rof

  } // rof

  TAddFilter::Pointer add = TAddFilter::New( );
  add->SetInput1( img1 );
  add->SetInput2( img2 );

  TImageWriter::Pointer output_image_writer = TImageWriter::New( );
  output_image_writer->SetInput( add->GetOutput( ) );
  output_image_writer->SetFileName( output_image_file );
  try
  {
    output_image_writer->Update( );
  }
  catch( itk::ExceptionObject& err )
  {
    std::cerr << "Error caught: " << err << std::endl;
    return( 1 );

  } // yrt

  input_prs.close( );

  return( 0 );
}

// eof - $RCSfile$
